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Abstract Recently, Dammak and coworkers [T] proposed that the quantum 
statistics of vibrations in condensed systems at low temperature could be 
Q ■ simulated by running molecular dynamics simulations in the presence of a 

c/3 I colored noise with an appropriate power spectral density. In the present con- 

tribution, we show how this method can be implemented in a flexible manner 
and at a low computational cost by synthesizing the corresponding noise 'on 
the fly'. The proposed algorithm is tested for a simple harmonic chain as well 
as for a more realistic model of aluminium crystal. The energy and Debye- 
I Waller factor are shown to be in good agreement with those obtained from 

harmonic approximations based on the phonon spectrum of the systems. The 
limitations of the method associated with anharmonic effects are also briefly 
' O , discussed. Some perspectives for disordered materials and heat transfer are 

Ch ' considered. 
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1 Introduction 
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Standard constant-energy molecular dynamics (MD) simulations explore the 
, phase space of a system by following the classical equations of motion, there- 

' fore producing a classical micro canonical distribution function. Very early 

I it was realized that these equations could be modified in several different 
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manners to produce a canonical phase space distribution at fixed temper- 
ature. Among the possible ways to achieve this is the implementation of a 
"Langevin thermostat" [2] that couples each particle independently to a fic- 
titious thermal bath by introducing in the equation of motion a random force 
and a friction term related by the classical fluctuation-dissipation theorem. 
This is illustrated by the following equation for a single degree of freedom X 
of mass M submitted to an external force F{X): 



M— = ~M-iV + F{X) + v/27\/7 e{t) (1) 

where V ~ dX/dt is the velocity, and is a white noise verifying: 

{0{t)9{t')) ^kBTSit^t'). (2) 

While this type of thermostat (and various other more sophisticated versions 
developed by Hoover, Nose and others, see for example Ref. U) provides 
a convenient control of the temperature, it is intrinsically classical, and in 
particular the equipartition theorem is strictly verified. Any quadratic degree 
of freedom will on average correspond to an energy kBT/2. While this is a 
minor problem in liquid systems at room temperature, in which quantum 
effects are effectively rather small, the situation is quite different in solids. 
The Debye temperature, that signals the onset of quantum effects, is often of 
the order of a few hundred Kelvins, meaning that quantum effects are indeed 
important. The most prominent manifestation of this fact is the deviation 
from the Dulong and Petit law of heat capacity, which is manifest in most 
solids already at room temperature. 

Obviously quantum dynamics in solids is well handled by the usual har- 
monic theory of solids, possibly including anharmonic corrections. Using this 
theory however requires a full calculation of the normal mode spectrum, 
which may be cumbersome for large or aperiodic systems. Calculations of 
transport properties such as heat conductivity also represents a real chal- 
lenge, as well as static properties in strongly anharmonic systems. While 
static aspects can be dealt with using path integral methods, dynamical ef- 
fects can not. Moreover, path integrals represent at low temperatures a rather 
costly alternative to much simpler classical simulations due to the large num- 
ber of replicas to include in the path integrals in order to converge to the 
quantum limit [l]. 

It is therefore highly desirable to develop methods that would allow one 
to describe correctly quantum effects at a relatively low computational cost, 
even if the representation is approximate. Recently two very interesting steps 
in this direction were taken by Dammak et al [T] and by Ceriotti et al [S]. 
The fundamental idea introduced in both papers is that the correct distri- 
bution of positions and impulsions for a quantum harmonic oscillator can be 
reproduced by using an equation similar to Eq. [1] but with a colored noise 
replacing the white noise d{t). In Ref. [S], the friction term was moreover 
substituted by a Non-Markovian term, thus permitting a greater flexibil- 
ity in the types of thermostats used. This non-Markovian term, and the 
corresponding noise, were replaced by a coupling to a finite set of classical 
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degrees of freedom forming the " quantum bath" , and the parameters charac- 
terizing the evolution of the quantum bath were fitted to obtain the correct 
energy distribution for a quantum harmonic oscillator. In Ref. [1], a slightly 
different approach was used. The friction term 7 was maintained Marko- 
vian, and the noise was replaced by a colored noise with power spectrum 
0{ijj) = fioj + {exp{fiuj/kBT) — In this manner, the correct energy 

distribution for the harmonic oscillator was directly obtained, and a very 
satisfactory agreement with low-temperature experimental data was demon- 
strated. The only drawback of the approach proposed in Ref. [1] is practical 
and concerns the colored noise generation. The authors used a method which 
is classical in the physics community [HIIZIIH]: but which requires to perform 
an inverse Fourier transform of the power spectrum of the noise over the 
entire time duration of the trajectory. The complete noise signal must there- 
fore be generated and stored in the computer memory before starting the 
simulation. The corresponding memory requirement, for large numbers of 
particles or for long trajectories, is important, making the implementation of 
the method rather cumbersome. Our contribution in the present paper is to 
show how this problem can be circumvented by using a synthesis of the noise 
employed in the signal processing community [HUH] , based on an appropriate 
linear filter applied to a white noise, thus making the quantum thermal bath 
method much easier to implement in practice. 

The manuscript is organized as follows. Section [2] recalls the basis of the 
method and discusses some physical issues related to convergence. Section [3] 
describes the synthesis of the noise using a linear filter. Applications to two 
selected examples are considered in Section |4l 



2 The quantum thermal bath method 

The principle of the quantum thermal bath is readily illustrated on the simple 
case of an harmonic oscillator. Using a unit mass and frequency loq, the 
equations of motion are 

X = V (3) 

V = -ujIx ~-fV + ^ e{t) (4) 

where X and V are the position and velocity of the oscillator, 7 the friction 
and 6{t) the external noise. The Fourier components arc easily obtained as 

X{^) = , . (5) 



and 



V{u:) = i^— . 6 



And finally the average energy of the oscillator is given by: 



(7) 
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where 0{u)) is the power spectral density (PSD) of the noise. In the foUowing, 
we wiU use capital letters and a "to denote spectral densities associated with 
a given random noise. 

For a constant 0{ijj) (white noise), integration of Eq. [7] shows that the 
energy is constant and independent of 7. If moreover 0(uj) = fcsT as required 
by the classical fluctuation-dissipation theorem, the energy is just equal to 
the thermal energy fcsT. 

For a colored noise with a finite frequency support [^Omax^ ^max], Eq. [7] 
yields, provided f2max > f^o+Ti E ~ 0{loo) with an accuracy that depends on 
7 (typically the actual value of the energy will average 0(u)) in the frequency 
range [cjq — 7/2,^0 + 7/2]). Hence, by choosing 

0(0.) = h\u:\(^+ ] ) (8) 

V2 exp(7i|w|/A;BJ ) — 1/ 

and a value of 7 small compared to ojq, one ensures that the energy of the 
oscillator is given by the Bose-Einstein formula including the zero point en- 
ergy. The reasoning is readily extended to a set of linearly coupled harmonic 
oscillators (since the friction and covariance of the noise are diagonal matri- 
ces). In this case, each eigenmode will have an energy equal to that of an 
harmonic oscillator with the corresponding frequency. 

A slight word of caution is in order here. The spectral density given in 
Eq. [5]does not have a finite support, and moreover it diverges proportionally 
to w for high frequencies. The integrand in Eq. [7] thus decays as l/ui for 
large w and the integral diverges at high frequency for any finite value of 7. 
Clearly this divergence is associated with the fact that the quantum thermal 
bath includes fiuctuations at arbitrary high frequencies, which are picked up 
by the broadening of the oscillator response function due to the damping 7. 
These high frequency fluctuations, however, are not expected to be present 
in any real system: in a crystal, the actual bath of oscillators that a given 
oscillator couples with is physically limited to frequencies of the order of 
the Debye frequency. No higher frequency is present in the system. Hence, it 
appears reasonable to modify Eq.|8]by introducing an upper cutoff frequency 
^max of the order of a few times the highest frequencies observed in the 
system. We will see in the following that such a cutoff frequency appears most 
naturally in a numerical implementation of the colored noise. The divergent 
contribution to the energy in Eq. [7]then behaves as jhi{Qrnax)-, and can be 
made arbitrarily small by choosing a small enough value of 7, once fimax has 
been specified. In general, and for a given cutoff frequency, the independence 
of the results on the damping parameter 7 should be checked to ensure that 
the relevant spectrum of frequencies is properly thermostated. 

On the other hand, it is also clear from Eq. [7] that, for a fixed cutoff 
^max, the result will be a decreasing function of the damping 7. In practice, 
7 should be kept small enough (compared to wq) that the integrand actually 
behaves like a S function, while the equilibration time remains reasonable. 
We will see in the following that values of j/loq ^ 10~^ provide a good 
compromise between these various requirements. 
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3 Noise synthesis 

The numerical implementation of the quantum bath as described above re- 
duces to generating a colored noise with a prescribed PSD, 0{uj). We use 
here a classical method of signal processing based on filtering a white noise 
signal j9j. For a continuous noise, we introduce the filter 

H{u:) ^ \f§{^) (9) 

and its inverse Fourier transform H{t) [O is given in Eq. The noise 0{t) 
is then obtained by convoluting H{t) with a random white noise r{t) with 
PSD R{uj) = 1. We have: 

/OO 
H{s)r{t - s)ds (10) 

The spectral density of the resulting noise is then \H [uj)]"^ R{uj) = 0{oj), 
which is the target PSD. This scheme can be implemented in a discrete al- 
gorithm suitable for computer calculation. The filter H{uj) is first discretized 
in 2Nf values with steps 5ui over an interval [—fimaxi ^max[- 

Hk = HikSuj), k = -Nf...Nf - 1. (11) 

The cutoff frequency fimax can be chosen on physical grounds as discussed 
above, and will typically be taken to be a few times the maximum eigenmode 
frequency of the system. A discrete Fourier transform yields H{t) over the 
interval [— 7r/(5w, 7r/(5w[ with a timestep h = T:/f2„iax- H{uj) being even and 
real, so is H{t): 

Hn = 7rTr Hk COS i^kn). (12) 

2Nf ^ ^Nf ^ ^ ' 

J k=-Nf 

The convolution is then performed discretely in the range [—tt/Suj, tt/Suj[, 
yielding a discrete colored noise with a time step h: 

Nf-l 

0n = 9{nh) = ^ HmTn-ra, (13) 
m=-Nj 

where is a white noise drawn from a gaussian or a uniform distribution 
with a variance \/h in order to yield a unity PSD. 

In practice, {Hm\m=-Nf...Nf-i is first computed using Egs.lTTIandfT^and 
stored. The noise generator is then initialized by drawing 2Nf initial values 
{T'm\m=o...2Nj-i and time is initialized to Nfh. Then, at each subsequent 
time step, with time nh and n > Nf, Vn-Nf is discarded and a new random 
number Tn+Nf is drawn and stored. The correlated noise is then obtained 
by computing the convolution in Eq. [T31 The present algorithm has a CPU 
computational cost close to a classical Langevin thermostat since the only 
difference is the convolution in Eq. 1131 The main difference is the memory 
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storage because the present algorithm requires to store ANf values per degree 
of freedom, corresponding to the filter and white noise. This requirement can 
however be easily managed and is much less than in the method used by 
Dammak et al [Tj , where an inverse Fourier transform is performed to create 
a signal over the complete duration of the simulation with a discretization 
that corresponds to the time step used by the MD integrator, i.e. a required 
storage (number of time steps) x (number of degrees of freedom). 

In principle, there is no specific reason to set the time step of the MD inte- 
grator (6t) equal to the time step of the noise generator h (which corresponds 
to the inverse of the cutoff frequency Umax)- Taking again the harmonic os- 
cillator as an example, a reasonable time step for integration of the equation 
of motion is St = O.Olo;,^^. The cutoff frequency for the noise, on the other 
hand, should be physically limited to a few times coq, let us say for concrete- 
ness f2max = 2a;o. Hence the noise time step can be taken to be M times 
larger than the MD time step, with M = 50 for this specific example. The 
simplest way to manage the two time steps is to generate the noise according 
to the procedure outlined above, with a time step h = M5t, and to keep 
the noise value constant for M integration time steps St. This amounts to 
generating a noise with an actual power spectrum 

S{n) = H{ojfC{oj)\ (14) 

where 

is the Fourier transform of the square function of amplitude unity and dura- 
tion h. Obviously the function S'(a;) is not the desired PSD 6?(w), but this 
problem is easily cured by replacing H[uj) in Eq. [Tl]by a corrected spectrum 

Hi{oj)^ H{uj)C{u:)-^ (16) 

so that the proper spectral density is generated over the interval [0, fimax]- 
As an illustration. Fig. [1] compares on a specific example the PSD of the 
noise generated by the above procedure with or without the final correction. 
Before correction, the intensity of the noise is below & and decreases rapidly 
at high frequencies. The reason is the suppression of the high frequencies 
due to the noise plateaus between successive updates. On the other hand, 
after correction, the PSD of the noise matches very well up to the cutoff 
frequency fimax- The results of the quantum thermal bath are illustrated 
below for two simple examples. 



4 Examples 

4.1 Linear chain of oscillators 

We consider first the very simple example of a linear chain of oscillators cou- 
pled by Hookcan springs between nearest neighbors, with periodic boundary 
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Fig. 1 (Color online) Noise generation at difTerent temperatures. The dashed lines 
show the target PSD at T = 10 K and T = 1000 K. The blue curves are the PSDs 
of the noise generated without correction and the red curves are after correction. 
The parameters of the quantum thermal bath used here are: fimax = 100 THz, 
Nj = 100, M = 30. 



conditions. The frequency of an individual spring is wq. In Fig. [2] the energy 
obtained from a MD run using a thermostat with flmax = 2a;o and 7 = ajo/50 
is compared to the exact result for this system. 

In Fig. [31 we show the influence of the friction parameter 7. As expected, 
a too large 7 results in a broadening of the oscillator response function which 
docs not pick up the correct range of frequency spectrum. As the spectrum 
has been cut off to a rather low frequency here, the result is an artificial 
decrease of the energy. However, a value of 7 of the order of wq/SO provides 
both a good approximation to the spectrum and reasonable equilibration 
time. 



4.2 Aluminium crystals 

We model an aluminum perfect crystal using the many-body embedded atom 
method potential developed by Ercolcssi and Adams [TU] . The simulation cell 
contains a perfect periodic crystal made of 6 x 6 x 6 repeating cells with 
crystallographic orientation X = [100], Y [010] and Z = [001]. The cell 
comprises N=864 atoms. We compute the evolution of the crystal energy and 
Debye- Waller factor as a function of temperature using the quantum thermal 
bath thermostat. In both cases, the results of the MD simulations are com- 
pared to classical mechanics calculations performed using a classical Langevin 
thermostat, to path-integral molecular dynamics (PIMD) calculations (see 
Ref [4] for details) and to quantum and classical harmonic approximations 
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Fig. 2 (Color online) Reduced energy (in units of hujo) as a function of reduced 
temperature Tiuo/ks; for a chain of 50 oscillators. The dots are the result from a 
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Fig. 3 (Color online) Reduced energy (in units of Tiuo) for reduced temperatures 
huo/kB ~ 0.2 and hujo/kB = 1.0 for a single oscillator, as a function of the friction 
coefficient ^/ujo- The point at 7 = is the exact energy. 
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Fig. 4 (Color online) Average energy per atom in an Aluminum crystal as a 
function of temperature. Different frictions 7 are used, averages over 5 x 10* MD 
steps after 5 x 10'' steps of equilibration. Parameters of the quantum thermal bath: 
Wmaa; = 101.34 THz, = 100, (5t = 1 fs, = 315t. Results are compared to quan- 
tum and classical harmonic approximations computed from the phonon spectrum, 
as well as PIMD simulations with a number P of replica satisfying PT = 4096 K. 



based on the phonon frequency distribution of the crystal, {wi}i=o;3Af-i- The 
latter is computed after diagonahzation of the Hessian matrix of the crystaL 
In the harmonic approximation, the quantum average energy is given by: 

3Ar-4 

(Eh 

4=0 

where the 3 translational modes have been omitted and O is the harmonic 
oscillator energy given in Eq. [8] The classical harmonic approximation is 
obtained by replacing O by fcsT. The comparison is shown in Fig. 01 The 
quantum thermal bath allows to recover very well the quantum harmonic ap- 
proximation at both frictions 7 = 0.1 THz and 1 THz and the figure shows 
the convergence to the classical regime and the gradual slight departure from 
the harmonic approximation at high temperatures. The highest frequency in 
the crystal being lOmax =59.7 THz, we conclude again that a ratio ^/iOmax 
in the range 1/50 to 1/500 is sufficient to ensure the convergence of the 
energy. In this example, Nf ~ 100, meaning that the memory requirement 
is about 100 times larger than would be the case for a classical molecular 
dynamics simulation, independently of the duration of the simulation. The 
present algorithm is therefore much less demanding than the previously pro- 
posed implementation [l] , in which memory storage was proportional to the 
duration of the simulation. 
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Fig. 5 (Color online) Debye- Waller factor, or average mean square displacement, 
in an Aluminum crystal as a function of temperature. The same parameters are 
used here as in Fig. [l] 



The second example is the Debye- Waller factor computed here as the 
atomic mean-square displacement (u^) = ;^ ~ ''i(O))^)- 1^. the 

harmonic approximation, the Debye- Waller factor is expressed as: 

, 3JV-4 ^ ^ 3Ar-4 A . ^ 

> - E ^ (^) ^ ^ E ikf 

i— 2—0 ^ 

The result is shown in Fig. [SJ where again a good agreement is seen be- 
tween the QTB MD simulations and the harmonic approximation at low 
temperatures. It should be noted however that, in contrast with the PIMD 
results, the QTB simulations tend to slightly overestimate the harmonic ap- 
proximation, a feature that will be discussed in section [5] below. 



5 Some words of caution 

In this section, we point out some differences between this method and the use 
of a classical Langevin thermostat. With a classical thermostat, the equations 
of motion can be transformed into a Fokker-Planck equation for the probabil- 
ity distribution of the positions and momenta. This Fokkcr-Planck equation 
admits as a stationary solution the classical Gibbs-Boltzmann distribution, 
independently of the manner in which each specific coordinate is coupled 
to the thermostat. In particular, with a set of weakly coupled (eg. trough 
a small anharmonic term in the energy) harmonic oscillators, the coupling 
of a finite number of oscillators to the thermostat is sufficient to ensure the 
thcrmalization of the whole system at the same temperature T. 
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Fig. 6 (Color online) Reduced energies (measured in units of huji) of two coupled 
quantum oscillators with frequencies oji and 0)2, as a function of the reduced tem- 
perature (measured in units of huji/kB)- The oscillator with the smaller frequency 
{(jji) is directly coupled to a thermal bath with the same parameters as used in 
figure [21 The second oscillator has a frequency UJ2 = V^lji and is not coupled to 
the thermal bath. The coupling between the two oscillators is ensured by a linear 
(squares) or quartic (dots) term He in the energy (see text). In both cases the 
strength of the coupling is e = 0.05 in reduced units. The points correspond to 
simulation data, the lines are the exact energies for decoupled oscillators. 



In the case of colored noise, the transformation to a Fokker Planck equa- 
tion is not possible, and there is to our knowledge no general expression 
available for the stationary phase space distribution. The quantum ther- 
mal bath presented above is coupled to all degrees of freedom, and ensures 
that each independent harmonic mode acquires the correct quantum energy, 
independently of the coupling between modes. One may however question 
the influence of such a coupling, and in particular the manner in which a 
" quantum" mode would transfer energy to other modes that are not directly 
coupled to the quantum bath. We have performed an exploratory numerical 
study of this situation by considering two weakly coupled oscillators, Xi and 
X2, with frequencies ioi and uj2 and identical masses, and coupled trough a 
weak coupling Hamiltonian Hc{Xi, X2). The first oscillator is coupled to the 
quantum thermal bath, while the second one evolves according to the classi- 
cal equations of motion, M X2 = — Mw|X2 — 8x2 He- We have considered the 
case of a linear coupling He — 6X1X2, and the case of a nonlinear coupling 
He = 6X^X2/2. In both cases, we find (see figure [5]) that both oscillators 
equilibrate, within numerical accuracy, at the desired energy. The equilibra- 
tion time in the case of a nonlinear coupling is however extremely long, and 
the error on the energy of the second oscillator remains significant even after 
10^ periods of oscillation. 
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In the more practical case of a weakly anharnionic system with a large 
number of modes, such as a crystal at low temperatures, a related ques- 
tion arises. Here each mode is coupled to the quantum bath with a coupling 
strength 7. However, the equation of motion is also affected by the coupling 
to other modes through the anharmonic terms of the energy. Classically, such 
anharmonic couplings give rise to a finite lifetime of the mode, characterized 
by an effective damping Fanh- In order for the quantum thermal bath to 
operate properly, it is desirable that the coupling strength 7 is stronger than 
this anharmonic coupling Fanh- On the other hand, as discussed previously, 
a large 7 also leads to inaccurate results for a single mode. One therefore has 
to find a reasonable compromise, in which the coupling is strong enough so 
that the mode can be thcrmalizcd within its lifetime F~^f^, but weak enough 
to keep the thermalisation accurate. In fact, the slight deviation observed for 
the Debye- Waller factor in Fig. [S] can be ascribed to the failure in finding a 
completely satisfactory compromise. By computing the lifetime of the short 
wavelength phonons using a classical simulation at a temperature that cor- 
responds to the quantum zero point energy, we find indeed that the order of 
magnitude of the lifetime is of a few picoseconds for the present interatomic 
potential. This implies that the friction 7 should be at least in the Terahertz 
range, which indeed yields the best results for the Debye- Waller factor. It is 
found that smaller values of 7 yield a repartition of energy among the normal 
modes of the crystal that tends to a classical distribution, independent of the 
mode frequency. 



6 Conclusions and perspectives 

In this manuscript, we have presented an "on the fly" implementation of 
a quantum thermal bath for molecular dynamics simulations, suitable for 
long calculations in large systems. Two simple test cases were considered 
as validations. The present algorithm is easily portable with limited memory 
requirements, and therefore opens the possibility of studying complex systems 
at a relatively modest computational cost. 

Direct applications of interest in this context include, for static properties, 
all cases in which the use of the harmonic approximation requires a numer- 
ically demanding matrix diagonalisation. This is the case in particular of 
amorphous or nanocrystalline systems at low temperatures. Using the quan- 
tum thermal bath allows one to access directly thermodynamic quantities 
from a simple molecular dynamics calculation. 

A more involved potential application, which may require further develop- 
ments, is the investigation of thermal transport properties in nanostructures 
and nanostructurcd materials, for which the use of classical approximations 
based on phonon spectra and lifetimes is not practical. The study presented in 
section [5] suggests that modes that are only indirectly coupled to the thermal 
bath also equilibrate correctly at their quantum energy. Therefore it should 
be possible, as is done in classical heat transfer simulations, to equilibrate a 
temperature profile by connecting a sample to two quantum thermal baths 
operating at different temperatures and in different regions of space. 
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Finally, the fundamental properties of the quantum thermal bath, and 
in particular of the phase space distribution created by the colored noise in 
a system which is not strictly harmonic, should also be investigated in the 
future. 
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